/*
 * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
 *
 * Please refer to the NVIDIA end user license agreement (EULA) associated
 * with this source code for terms and conditions that govern your use of
 * this software. Any use, reproduction, disclosure, or distribution of
 * this software and related documentation outside the terms of the EULA
 * is strictly prohibited.
 *
 */



#include <math.h>



///////////////////////////////////////////////////////////////////////////////
// Polynomial approximation of cumulative normal distribution function
///////////////////////////////////////////////////////////////////////////////
double cnd(double d){
    const double       A1 = 0.31938153;
    const double       A2 = -0.356563782;
    const double       A3 = 1.781477937;
    const double       A4 = -1.821255978;
    const double       A5 = 1.330274429;
    const double RSQRT2PI = 0.39894228040143267793994605993438;

    double
        K = 1.0 / (1.0 + 0.2316419 * fabs(d));

    double
        cnd = RSQRT2PI * exp(- 0.5 * d * d) * 
        (K * (A1 + K * (A2 + K * (A3 + K * (A4 + K * A5)))));

    if(d > 0)
        cnd = 1.0 - cnd;

    return cnd;
}


////////////////////////////////////////////////////////////////////////////////
// Process an array of optN options
////////////////////////////////////////////////////////////////////////////////

// extern "C" 

void BlackScholesLLC(
                float *h_CallResult,
                float *h_PutResult,
                float *h_StockPrice,
                float *h_OptionStrike,
                float *h_OptionYears,
                float Riskfree,
                float Volatility,
                int optN
                ){
        #pragma acc kernels loop name("blackscholes") copyin(h_StockPrice, h_OptionStrike, h_OptionYears, Riskfree, Volatility) copyout(h_CallResult, h_PutResult)
        for(int opt = 0; opt < optN; opt++)
        {
                // Black-Scholes formula for both call and put
                double S = h_StockPrice[opt], X = h_OptionStrike[opt], T = h_OptionYears[opt], R = Riskfree, V = Volatility;
                double sqrtT = sqrt(T);
                double    d1 = (log(S / X) + (R + 0.5 * V * V) * T) / (V * sqrtT);
                double    d2 = d1 - V * sqrtT;
                double CNDD1; // = CND(d1);
                double CNDD2; //  = CND(d2);
                //Calculate Call and Put simultaneously
                double expRT = exp(- R * T);

                {
                        CNDD1 = cnd(d1);
                        CNDD2 = cnd(d2);
                        h_CallResult[opt]   = (float)(S * CNDD1 - X * expRT * CNDD2);
                        h_PutResult[opt]   = (float)(X * expRT * (1.0 - CNDD2) - S * (1.0 - CNDD1));
                }

        }
}
